Chaos in Pseudo-Newtonian Black Holes with Halos. 
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The Newtonian as well as the special relativistic dynamics are used to study the stability of 
orbits of a test particle moving around a black hole plus a dipolar halo. The black hole is modeled 
by either the usual monopole potential or the Paczyhki-Wiita pseudo-Newtonian potential. The full 
general relativistic similar case is also considered. The Poincare section method and the Lyapunov 
characteristic exponents show that the orbits for the pseudo-Newtonian potential models are more 
unstable than the corresponding general relativistic geodesies. 
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I. INTRODUCTION. 

o 

To consider relativistic effects in many body simulations is not a simple task due to the fact that the metric 
representing their gravitational interaction is far from being known. For the simplest case of two gravitating bodies 
the metric is known numerically only for few initial conditions and for a limited amount of time [see for instance, 
Marronetti et al. ( [|2000f| . Also, assuming that the metric is known, the use of the geodesic equations to determine 
' 4 ' the trajectory of the bodies represents a quite non trivial problem. 

0^ \ In general, we have three main ways to consider complex systems: a) A full numeric approach with its inherent 
limitations due to the use of floating point arithmetics and arbitrariness of discretizations of fundamentally continuous 
functions and variables. Also we have rather unphysical ad hoc assumptions like the introduction of numerical viscosity, 
b) The use of perturbative methods that are usually employed together with drastic approximations like the mean 
field approximation for the potentials in many body simulations. These approximations introduce irreversibility in an 
_^ , intrinsic reversible situation, c) The modeling of the problem with simpler equations in which one takes into account 
-H ' a few essential features of the problem. In general, this model can be solved in a more exact form of the two precedent 
case. But, we have changed the initial problem for a simpler one that may falsify results. In other words, there is not 
• a perfect method to solve a a complex problem. We believe that all of them are valid when the adequate cautions are 
taken. Furthermore, they are complementary and the usually not proven mathematical or internal consistence of the 
Q . methods can be independently checked at least for some particular cases. 

Oh Due to the weakness of the gravitational field, far from the particles' horizon, the Newtonian gravity is proven to 
be a reliable description of the gravitational interaction. One can simulate relativistic effects within the Newtonian 
theory changing the usual potentials to take into account the existence of the horizon. In other words to model 
55 \ relativistic effects via a pseudo-Newtonian potential. These models are simpler enough to describe complex systems 
that are far beyond of todays knowledge of the full general relativity, e.g., the n-body simulation of the collision of 
two galaxies to any degree of resolution. 

One of the simplest pseudo- Newtonian potentials to describe behavior of test particles moving close to a black hole 



X 



is the Paczynski and Wiita ( [ 198C ]) pseudo-potential 



* = "#^. (1) 

R Rg 

The addition of the term R g = 2GM/c 2 critically changes the particles' trajectory near the source. Some results like 
the last stable circular orbit are predicted in this model. Other pseudo-Newtonian models can be found in literature, 



e.g., the one studied by Semerak and Karas ( [ 1999 1 ) to describe rotating black holes, i.e., to approximate the Kerr 
solution. 

We believe that the study of the Paczynski and Wiita (PW) potential in simple albeit nontrivial situations may shed 
some light into the correctness of the pseudo-potential approach. In particular, in this article, we study integrability 
and chaos in a system that represents a spherically symmetric source (monopole) surrounded by a dipolar halo (external 
dipole), that is the simplest mean potential used to describe astrophysical systems restricted to a core and halo, see 



for instance Binney and Tremaine ( [1987]). Different theoretical approaches are used to study this configuration. 
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First we use Newton second law to find the motion equations for test particles (a = — V$) for two different potentials 
that describe a core plus a dipolar halo system: a) The standard monopole plus external dipole expansion that solves 
the usual Laplace equation that is totally integrable, see for instance Grammaticos et al.( | 1985[ |), and b) We replace 
in the former case the monopole term by the PW potential (|l|) . In this case the trajectories are chaotic like in the 
equivalent full general relativistic system, Vieira and Letelier ( ]1997|] ). 

We also analyze the equivalent cases using the special relativistic dynamics. We solve the equation — with 
a' 1 = = ljf t (t^jt) an d = j(— V$-v/c, — V<£>), where 7 = (I — v 2 /c 2 ) -1 / 2 , and $ is taken as in the Newtonian 
cases. We first use the monopole plus dipole potential that solves the Laplace equation. A phase space analysis shows 
that the system is stable. Replacing the monopole term by the PW potential we obtain a very unstable system. 
We also review the equivalent system in general relativity. The geodesic equations for Schw arzsch ild monopole plus 
dipolar halo give us chaotic trajectories in the phase space as shown in Vieira and Letelier ( [ 1997 1 ) . 

In each one of the studied cases we have an integrable Hamiltonian system of equations for the motion of a test 
particle moving in a spherically symmetric attraction center (standard monopole, PW potential or Schwarzschild 
metric) that is perturbed by an external dipole ter m. In all these situations we can apply the KAM (Kolmogorov, 
Arnold and Moser) theory, see for instance Tabor ( [1989]). Since our mass distribution has axial symmetry we are 
restricted to an effective two-dimensional problem. In the integrable case, in phase space, the orbits of test particles 
will be confined to a 2-torus. For a constant value of one of the coordinates we obtain a planar section of the phase 
space. In the integrable case we shall see closed curves for each initial condition, intersections of invariant tori. While 
in the non-integrable case some tori will be destroyed and the region will be ergodically fulfilled. In order to evaluate 
the degree of instability of the orbits in each system we also compute the Lyapunov exponents that indicate us how 
initially close trajectories separate. 



II. NEWTONIAN DYNAMICS. 



The standard monopole plus external dipole potential in the usual cylindrical coordinates (r, z, <p) is 

GM 

$ = +Dz, (2) 

yr 2 + z 2 

where D is the dipolar strength, G the Newton constant, and M the mass of the attraction center. We use units such 
that GM — 1, furthermore we shall take c = 1. From the angular momentum and energy conservation we find that 
the motion is restricted to the region defined by 

L 2 

E 2 - 1 - -j - 2$ > 0. (3) 

L is the specific angular momentum of the test particle and E = vT + 2E mec h, where 

f 2 -+- 7 2 L 2 

E m ech = + *(r, Z) + ^ 

is the specific energy. Note that E become imaginary for Emech < —0.5 that is the energy of a particle standing on 
the black hole horizon. The phase space orbits are studied using the Poincare section method. In Fig.l we present the 
surface of section z = for the constants: L = 3.9, E = 0.976, and D = 2 x 10~ 4 . This surface section characterizes 
an integrable system as expected. 

Now we shall replace the monopolar term by the PW pseudo-Newtonian potential, i.e., 

$ = 1 +Dz. (4) 

Vr 2 + z 2 -2 

Again, the motion of test particles will be restricted to the region that solves (||) with $ given by (||). In Fig. 2 we 
present the surface of section z = 0. We take the values for the constants as in the preceding case: L = 3.9, E — 0.976, 
and D = 2 x 10~ 4 . Contrary to the previous case we observe chaotic orbits in this Poincare section. 
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III. SPECIAL RELATIVISTIC DYNAMICS. 



In principle, the use of the special relativistic dynam ics sh ould improve the modeling of general relativity with 
pseudo-Newtonian potentials, see Abramowicz et al. ( [ 1996j ]). Although, these authors found that the predicted 



spectra often differ rather substantially from those obtained in the full general relativity context. ^From the relativistic 
motion equation we get 

^(7 + *) -0^7 + *-^, (5) 
dO _ L 
dt 7 r 2 

By using the above equations and u^u^ = 1 we obtain 

{E -$) 2 (l-f 2 - z 2 



(G) 



L 2 



[{E-<f>)r} 2 ' 

which is used to calculate the region in which the motion is confined. Finally, the motion equations for the variables 
r and z are, 



($- 






dr 2 
dt ' 


(*- 


d 2 z 




dz 2 
~dt ' 



d<& dz dr L 2 
~dz~dtdt ~ (E - $)r 3 
<9<I> dz dr 
dr dt dt 



(7) 



(8) 



As in the previous section, we start with the usual monopole plus external dipole potential field, i.e., we identify <!> 
with (0). In Fig. 3 we draw the Poincare section defined by the plane z = 0. The constants are the same of the 
preceding section, L = 3.9, E = 0.976, and D = 2 x 10~ 4 . We notice that the tori were preserved in this case, we 
have stability of orbits. This is and indication of integrability of the system. 

Now we start the study of the PW potential plus dipolar halo, i.e., we identify $ with Unfortunately we cannot 
confine the orbits by using the constants attributed to all the preceding cases. We put, L = 4.2, E = 0.972, and 
D = 4.2 x 10~ 4 . Now the Poincare section is taken as z = —5. The figure in this case, Fig. 4, represents a very 
chaotic system. We used the same constants to draw another Poincare section for PW potential plus dipolar halo 
using Newtonian dynamics. The results are presented in Fig. 5. We see some stable islands in the negative p r region 
that cannot be observed Fig. 4. We conclude then that the orbits obtained in the special relativistic context are less 
stable than the ones obtained with Newton law. The conjugated variables used were dr/dt and r. We made some 
tests using dr/dr and r. The results were qualitatively the same. 

IV. GENERAL RELATIVISTIC DYNAMICS. 

We start from the axisymmetric line element 

ds 2 = e W(u,v) dt 2 _ e -W(u,v)( u 2 _ _ v 2 )d(j) 2 (Q) 

_ e 2( 7 («,,)-^)) (u 2 _ v 2 } ( + 

\u 2 — 1 1 — V 2 J 

in prolate coordinates (t,u,v ,</>). The coordinates u and v are related to the usual cylindrical coordinated by u — 
(R + + i?_)/(2m) and v = (R + - R-)/(2m), where R± = [r 2 + (z ± m) 2 ] 1 / 2 and m= GM/c 2 . The Schwarzschild 
monopole plus a dipolar halo is represented by 

ip(u,v) = ilog (y - ^) +Duv. (10) 

Note that taken the limit, lim c -2 =0 ip/c" 2 , with aid of PHopital rule, we recover (||). To have the right units to take 
the limit we need to add a c~ 2 factor to D. 

The Einstein equations for these class of solutions as well as the corresponding geodesic equations are studied in 
great detail in Vieira and Letelier ( |199£ |). Due to the axial symmetry of the metric again the effective geodesic 



dynamics of the test particles is restricted to a three dimensional "phase space" 
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The Poincare section is draw for v — (that is equivalent to z = 0). In Fig. 6 we present the section for the values 
of the constants L = 3.9, E = 0.976, and D = 2 x 1CT 4 . Chaotic orbits may be observed for instance in the region 
indicated with a rectangle. A zoom of this region is presented in Fig. 7. We can compare Fig. 6 with Fig. 2 and 
conclude that the orbits obtained via geodesic equation in general relativity are more stable than the ones obtained 
from the PW potential plus dipolar halo in Newtonian and special relativistic dynamics. 



V. LYAPUNOV EXPONENT 



We shall study the Lyapunov exponents for the systems above described to better analyze the orbits stability. We 
shall use the Lyapunov characteristic number (LCN) that is defined as the double limit 



LCN = lim 

<5o -» 



log(V<5 ) 



(11) 



where 5q and 5 are the deviation of two near by or bits at times and t respectively. We get the largest LCN using 



the technique suggested by Benettin et al. ( [1976] ) 

We start comparing the LCN for orbits in a PW+Dipole system in special relativity and the LCN for orbits in a 
PW+Dipole in Newtonian theory. The Constants are L = 4.1, E = 0.972, and D = 4.1 X 10~ 4 . The maximum LCN 
was obtained around r = 20, z = —5, and p r = —0.04. Note that the value of p z is determined by the constants of 
motion and the value of r, z and p r . For the relativistic case we get LCN = (3.2 ± 0.4) x 10~ 4 while for the Newtonian 
approach we obtain LCN = (1.8 ± 0.4) x 10~ 4 . We did some tests for the usual integrable Newtonian monopole plus 
dipole system and we always obtain for the LCN at least one order of magnitude lower that the precedent cases. 

For orbits of test particles in the the full general relativistic monopole plus dipole system and in the Newtonian 
PW+Dipole system we chose L = 3.902, E = 0.9756 and D = 2.0 x 10~ 4 . We obtain for orbits in the PW+Dipole 
system LCN = (2.0 ± 0.5) x 10~ 4 . This value was obtained for orbits around r = 7.5, z = 0, and p r = 0. For the 
general relativistic system the proper time and the coordinate time were tested in the equation ( p+| ) and no significant 
difference was found. The largest LCN was computed around u = 9.75, v = 0, &ndp u = —0.038. As before, p v is fixed 
by the value of the other variables and the motion constants. We found always LCN < 5 x 10 5 . The Lyapunov like 
coefficients used in general relativistic systems may have different forms as the one suggested by Burd and Tavakol ( 



[1993 1) in the study of Bianchi IX systems. However, we have studied a simple system with no singularities besides the 
black hole where we have a well defined evolution parameter. Hence, in this case no significant difference should be 
found by using other definition of the Lyapunov coefficients. Furthermore, in the general relativistic system studied 
we have several natural ways to choose the space variables e.g., the spheroidal {u, v, ip) an d the cylindrical (r, z, <ft). We 
found no significant differences when either system of coordinates are used to describe the orbits of particles moving 
a few Schwarzschild radii apart from the central black hole. 

In summary, the study of Lyapunov coefficients confirms the qualitative analysis of the Poincare section method, 
we have that the general relativistic orbits are more stable than the Newtonian and special relativistic ones. The 
special relativistic orbits are the most unstable. 



VI. DISCUSSION 



In the Paczyhski-Witta potential, the term —2GM/c 2 in the denominator of the equation (|J) creates a saddle point 
in the effective potential in Newtonian as well as in special relativistic dynamics. The addition of the dipole term 
se parat es the stable and unstable manifold emanating from the hyperbolic fix point as discussed by Letelier and Vieira 
( [ 1998 [). I n this case, as consequence of the Poincare-Birkhoff theorem, th ere is an homoclinic web that gives rise to 
chaotic motion for bounded orbits in phase space, see for instance Tabor ( |198£]). 

The chaotic orbits encountered in the pseudo-Newtonian plus dipole system agrees with the general relativistic 
equivalent situation. However, those effects might be distorted in the PW approach because the Poincare sections as 
well as the Lyapunov exponents show more unsta ble orbits. This instability is magnified when the special relativistic 
dynamics is used. Vokrouhlicky and Karas [1998] studied the stability of orbits for parti cles gravitating around a 1/ R 
Newtonian potential with an axisymmetric perturbation. Sridhar and Touma ( ] 1999| ) found for the same class of 
potentials that the instability decreases in orbits closer to the black hole. This result may not be verified when pseudo 
Newtonian or full general relativistic model are considered. The main difference being the presence of a saddle point 
in the effective potential near the black hole. Therefore orbits near the core may be more unstable because of this 
critical point in the effective potential that is a source of instability. 



4 



ACKNOWLEDGMENTS 

The authors thank to FAPESP for financial support. EG thanks A.E. Motter for discussions about Poincare-Birkhoff 
theorem. 



1996 
1976 
1987 
1993 
1985 
1998 
2000 
1980 
1999 



Abramowicz M.A., Beloborodov A.M. , Chen X.M., Igumenshchev I.V., 1996, A&A 313, 334. 
Benettin C, Galgani L. and Giorgilli A. , 1976, Phys. Rev. A14, 2338. 
Binney J., Tremaine S., 1987, Galactic Dynamics, Princeton University Press. 
Burd A and Tavakol R, 1993 Phys Rev. D 47, 5336. 

Grammaticos B., Dorizzi B., Ramani A., Hietarinta J., 1985, Phys. Lett. A 109, 81. 
Letelier P.S., Vieira W.M., 1998, Phys. Lett. A 242, 7. 

Marronetti P., Huq M., Laguna P., Matzner R.A., Shoemaker D., 2000, Phys. Rev. D 62, 401. 
Paczyhski B., Wiita P.J. ,1980, A&A 88 ,23. 
Semerak O., Karas V. , 1999 , A&A 343 , 325. 
1999 Sridha S. and Touma J., 1999, MNRAS 303, 483. 

Tabor, M., 1989, Chaos and Integrability in Nonlinear Dynamics, John Wiley&Sons , New York. 
Vokrouhlicky D. and Karas V., 1998, MNRAS 298, 53. 
Vieira W.M., Letelier, P.S., 1997, Phys. Lett. A 228, 22. 
Vieira W.M., Letelier, P.S., 1999, ApJ., 513, 383. 



1998 
1997 
1999 



5 



FIG. 1. Surface of section for the Newtonian motion of 
a test particle in a standard monopole plus external dipole 
potential for L z = 3.9, E = 0.976, and D = 2 x 1CT 4 . The 
section corresponds to the plane 2 = 0. For these values of 
the parameters we have the section of an integrable motion. 



FIG. 4. Surface of section for the special relativistic mo- 
tion of a test particle in a Paczyriski-Wiita potential plus a 
dipolar halo for L z = 4.2, E = 0.972, and D = 4.2 x 10~ 4 . 
The section corresponds to the plane z = — 5. We have a very 
irregular motion. 





FIG. 2. Surface of section for the Newtonian motion of a 
test particle in a Paczyhski-Wiita potential plus a dipolar 
halo for L z = 3.9, E = 0.976, and D = 2 x 10~ 4 . The section 
corresponds to the plane 2 = 0. We see chaotic motion. 




FIG. 3. Surface of section for the special relativistic motion 
of a test particle in a usual monopole potential plus a dipolar 
halo for L z = 3.9, E = 0.976, and D = 2 x 10~ 4 . The 
section corresponds to the plane 2 = 0. For these values of 
the parameters we have the section of a regular motion 



FIG. 5. Surface of section of the Newtonian motion of a 
test particle in a Paczyriski Wiita potential plus a dipolar halo 
for L z = 4.2, E = 0.972, and D = 4.2 x 10~ 4 . The section 
corresponds to the plane z = —5. We have an irregular motion 
but it is more stable than the one shown in the precedent 
figure. 




FIG. 6. Surface of section for the geodesic motion of a test 
particle in a Schwarzschild monopole with a dipolar halo for 
L z = 3.9, E = 0.976, and D = 2 x 10~ 4 . The section cor- 
responds to the plane v = 0. For these parameters we have 
small regions of instability. 



FIG. 7. A zoom of the small boxed region of the previous 
figure. 
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